clc; clear; close all;

ugrid = [0.01, 0.1];

for u = 1:2
    
upsilon = ugrid(u);

sigma = 4; theta = 3.8;  alpha = 1;

I = 3; K = 3000;
 
%%%%%%%%%%%%%%%%%% Trade Fundamentals %%%%%%%%%%%%%%%%%%%%%%
b = 1*ones(I,1);
tau =   ones(I,I);
f =  0.02.*ones(I,I);
for i = 1:I
    tau(i,i) = 1;
    f(i,i) = 0.02;
end
L = 10*ones(I,1);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
GG = 70;
grid = exp(linspace(log(9),log(1),GG));

for g = 1:GG 
    
tau = grid(g)*ones(I,I);
for i = 1:I
   tau(i,i) = 1; 
end

w = ones(I,1);
y = ones(I,1);
P = ones(I,1);
wr = ones(K,I);
rho = (linspace((10^(-5)),(1-10^(-5)),K))';
 
tol = 1; iter = 0;

while tol > 0.0001
    
iter = iter + 1;

V = (sigma./(sigma-1).*tau).^(1-sigma).*repmat(P',I,1).^(sigma-1).*repmat(L'.*y',I,1);

tolw = 1;

while tolw > 0.001

for i = 1:I
    for j = 1:I
    obj = abs((1-rho).^((1-sigma)./theta).*wr(:,i).^(1-sigma) - b(i).^(sigma-1).*sigma.*w(j).*f(i,j).*V(i,j).^(-1));
    
    star(i,j) = rho(find(obj==min(obj),1));
    pi{i,j} = max(b(i).^(sigma-1).*(1-rho).^((1-sigma)./theta)./sigma.*V(i,j).*wr(:,i).^(1-sigma) - w(j).*f(i,j),0);
    rev{i,j} = (pi{i,j} >= 0).*(b(i).^(sigma-1).*(1-rho).^((1-sigma)./theta).*V(i,j).*wr(:,i).^(1-sigma));
    end
end

PI = zeros(K,I); REV = zeros(K,I);
for i = 1:I
    for j = 1:I
    PI(:,i) = PI(:,i)+pi{i,j}; 
    REV(:,i) = REV(:,i) + rev{i,j};
    end
end

new_wr = repmat(w',K,1).^(1-upsilon).*PI.^(upsilon);
for i = 1:I
   new_wr(:,i) = new_wr(:,i)./trapz(rho,new_wr(:,i)); 
end

tolw = sum(sum(abs(new_wr-wr)));
wr = 0.1*new_wr + 0.9.*wr;

end

for i = 1:I
   d(i,:) = trapz(rho,PI(:,i)); 
end

yr = wr + PI;

for i = 1:I
   y(i,:) = trapz(rho,yr(:,i)); 
end

for i = 1:I
    for j = 1:I
    pp(i,j) = trapz(rho, ((1-rho).^((1-sigma)./theta).*wr(:,i).^(1-sigma)).*(pi{i,j}>0));
    end
end

X = repmat(L.*b.^(sigma-1),1,I).*V.*pp;
lambda = X./repmat(sum(X,1),I,1);
 
newP =  sum(repmat(L.*b.^(theta-1),1,I).*(sigma./(sigma-1).*tau).^(1-sigma).*pp,1).^(1./(1-sigma));
newP = newP';

tol = sum(abs(newP-P));
P = 0.9*P + 0.1*newP;


end

YY(:,g) = yr(:,1);


%=========================================================================%
Average = trapz(rho,yr(:,1)); 
Gini(g,:) = trapz(yr(:,1),(1-rho).*rho)./Average(1,:); 
Gini_D(g,:) = sum(sum(abs(repmat(wr(:,1),1,K)-repmat(wr(:,1)',K,1))))./(2.*K.^2.*mean(wr(:,1)));

rhostar(g,:) = 1-star(1,2);

onerho = (1-star);
Gini_pred(:,g) = alpha.*(sigma-1)./(sigma*theta) - 2.*alpha.*(theta-(sigma-1))./(sigma*theta).*(sigma-1)./(4*theta-2*(sigma-1))...
    .* sum(lambda.*onerho.*repmat((L.*y)',I,1)./repmat(L.*y,1,I),2);

end

Gini_pred = Gini_pred(1,:)';

RR(:,u) = rhostar;
GG_pred(:,u) = Gini_pred;
GG_actual(:,u) = Gini;

clearvars Gini_pred rhostar Gini

end

plot(RR(:,1), GG_actual(:,1)./GG_actual(1,1),'--', RR(:,1), GG_pred(:,1)./GG_pred(1,1),'-'); hold on;

plot(RR(:,2), GG_actual(:,2)./GG_actual(1,2),'--', RR(:,2), GG_pred(:,2)./GG_pred(1,2),'-'); 



 





 